Reconstructing quantum states efficiently 
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Quantum state tomography, the abiUty to deduce the density matrix of a quantum system from measured data, 
is of fundamental importance for the verification of present and future quantum devices. It has been reaUzed 
in systems with few components but for larger systems it becomes rapidly infeasible because the number of 
quantum measurements and computational resources required to process them grow exponentially in the system 
size. Here we show that we can gain an exponential advantage over direct state tomography for quantum states 
typically realized in nature. Based on singular value thresholding and matrix product state methods we introduce 
a state reconstruction scheme that relies only on a linear number of measurements. The computational resources 
for the postprocessing required to reconstruct the state with high fidelity from these measurements is polynomial 
in the system size. 



It is one of the principal features distinguishing classical 
from quantum many-body systems, that for the former the 
specification of a state requires a parameter set whose size 
scales linearly in the number of subsystems, while in the lat- 
ter this set scales exponentially. It is this difference that sup- 
ports the observation that quantum devices appear to be fun- 
damentally hard (exponential in the number of subsystems) to 
simulate on a classical computer and may in turn possess a 
computational power exceeding that of classical devices HI. 

This presents us with the blessing of being able to construct 
information processing devices fundamentally superior to any 
classical device and the curse of their complexity, challenging 
our ability to verify efficiently that such a quantum informa- 
tion processing device or quantum simulator is actually func- 
tioning as intended. Such devices and, more generally, quan- 
tum simulators may be viewed as the preparation of elaborate 
quantum states on which we then carry out measurements. 
Verifying efficiently that an intended state — the ground or 
thermal state of a physical Hamiltonian for example — has in- 
deed been prepared by a quantum information processor or a 
quantum simulator is essential. 

The full determination of the quantum state of a sys- 
tem, that is quantum state tomography |2|, can of course be 
achieved - one simply measures a complete set of observ- 
ables whose expectation values fully determine the quantum 
state IsHtI. In practice however, this approach is beset with 
several problems when applied to quantum-many party sys- 
tems. Firstly, in quantum state tomography the size of the set 
of measurements scales exponentially with the number of sub- 
systems. For moderately sized systems, such as the electronic 
state of 8 ions |4 |, tomography has been demonstrated but it 
rapidly becomes infeasible for larger systems thanks both to 
excessive time required to carry out the measurements and be- 
cause the precision of those measurements has to increase ex- 
ponentially to ensure that a function of the probability ampli- 
tudes of the state will not return an essentially random result. 
Secondly, making the connection between the measurement 
data on the one hand and the density matrix of a state best ap- 
proximating these data on the other will usually require clas- 
sical postprocessing that cannot be executed efficiently on a 
classical computer (see |4|). Thirdly, writing out the full state 
of a physical system will be impossible for more than approx- 
imately 40 spin- 1/2 particles and approximate representations 



from which one can extract expectation values efficiently to 
high precision need to be used. 

Here, we address all of the above challenges at the same 
time and demonstrate the efficiency of the proposed approach 
with numerical examples. We present the case of general pure 
quantum states in some detail and outline generalizations to 
mixed states. To pave the way towards the general argument, 
we start the exposition with a discussion of unique ground 
states of local Hamiltonians. 

Consider a /c-local Hamiltonian acting on a set LofN spins 
arranged on some lattice equipped with the notion of a dis- 
tance dist{i^ j) between sites i and j, 

H = Y^k, (1) 

where each Hamiltonian hi acts on spins that are at most 
at a distance k from spin i. Let us collect these in the set 
li = {j e L : dist(z,j) < k}. The non-degenerate ground 
state of a /c-local Hamiltonian is the state of lowest energy 
and therefore uniquely determined by the expectation values 
of the hi. Indeed, if there was another state with the same 
expectation values, its energy would be the same, violating 
the uniqueness assumption. For an unknown /c-local Hamilto- 
nian we do not know the hi and cannot restrict measurements 
to these observables only. The ground state \gs) is neverthe- 
less uniquely determined by all its reductions to the sites li, 
Qi = tr L\j.[\gs){gs\], as these determine all possible expec- 
tation values of operators acting on li, in particular those of 
the unknown hi. Hence, if an experiment prepares the unique 
ground state of some unknown /c-local Hamiltonian, we can 
determine that state fully by measuring the N reduced density 
matrices Qi. In this setting we have hence overcome the first 
problem mentioned above as the state is fully determined by 
the Qi - so N density matrices of size 2 1^*1 x 2 1^*1, where \Ii\ 
is the cardinality of li . For a nearest neighbour Hamiltonian 
on a d-dimensional lattice, e.g., |/^| = 2d -\- 1. Using these 
insights, we will discuss the case of general pure states later. 

To overcome the second problem, we require an efficient 
method to find a pure state {ip) whose reduced density matri- 
ces ai = tr coincide with the Qi. The method of 
choice is singular value thresholding (SVT) 1 8 , 9 1 (see the Ap- 
pendix for technical details), which has been developed very 
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recently in the context of classical compressive sampling or 
matrix completion and may also be applied to the quan- 
tum setting (lIlIISl. SVT provides a recursive algorithm that 
converges provably towards a low rank solution satisfying a 
set of linear constraints such as the requirement to match re- 
duced density matrices. SVT converges rapidly to the solution 
especially so when it has low rank. Unfortunately, SVT as 
originally proposed is not scalable as it requires a full repre- 
sentation of a matrix of the same size as the density matrix 
describing the system and a singular value decompositions 
of this matrix. Hence, both the requirement for memory and 
time scale exponentially in the number of sub- systems and the 
straightforward application of SVT is restricted to well below 
20 spin- 1/2 particles. However, as we will see, a modification 
of the algorithm allows us to overcome this problem. 

Denote by \(j)) the unknown target state and hy Qi, i = 
1, . . . , A^, its reduced density matrices as described above. 
The following modification of the standard SVT algorithm 
(see appendix) will yield a state {ip) whose reduced density 
matrices closely matches those of Let cr^, a = 

X, 0, the Pauli spin matrices acting on site j and by 
denote all the possible operators Ylj^j. cr^^ , i = 1, . . . , A^, of 
which there are XliGL^'^'' ~- Fi'om Qi we know the ex- 
pectation values (01 n^.^^^ a^' |(/)) = tr/. [gi Ujeh 

hence the numbers pk = ((/>|^fe|</>), ^ = 1, . . . ,iir. The al- 
gorithm may then be described as follows. First set up the 
operator R = Y^^^iPkPk/"^^ and initialize Yq (e.g., by the 
zero matrix). Then proceed inductively by finding the eigen- 
state \yn) with largest eigenvalue, i/n, of Yn and set 

/c=l 

(2) 

A rigorous proof of convergence of = \xi.[\yn){yn\] to 
Qi (equivalently of {yn\Pk\yn) to pk) will be presented else- 
where. Heuristically, convergence is suggested by the exten- 
sive numerics below and can be expected from the fact that 
SVT possesses a convergence proof for small (5r, G R 191, see 
Appendix. 

So far, this algorithm still suffers from the fact that in ev- 
ery step the 2^ x 2^ matrix Y^ needs to be diagonalized. 

However, the Yn are of the form Ylk=i ^kPk, CLk ^ R, i-^., 
they have the form of a local "Hamiltonian". In one spatial di- 
mension, the ground states of local Hamiltonians are well ap- 
proximated by matrix product states (MPS) |1^ J4J. Hence, 
\yn) can be determined employing MPS algorithms |[T5HT6ll , 
for which the number of parameters scale polynomially in the 
system size and converge rapidly LIZ. JM • Hence, for one 
spatial dimension, we have overcome the second and third 
problem mentioned above: This postprocessing is efficient as 
MPS algorithms are and a MPS provides an efficient represen- 
tation of the state as it depends only on linearly many parame- 
ters. Any general pure state may be represented by a MPS and 
generic MPS are unique ground states of local Hamiltonians 
ifTTI . Hence, the above algorithm will produce a state that is 
close to the target state if the MPS dimension and the size of 



the reduced density matrices is chosen sufficiently large. 

For higher spatial dimensions MPS are not efficient repre- 
sentations and for optimal performance they need to be re- 
placed by other variational classes. A variety of MPS gener- 
alization have been proposed of which the most promising are 
perhaps the tensor- tree ansatz [il9l and MERA-approach 1201 , 
PEPS \2V\ and weighted graph states |22|. For each of these, 
numerical algorithms have been developed that determine the 
largest eigenvalue of a local Hamiltonian. Being the key in- 
gredient in our modified SVT method developed here, these 
more general variational classes may be combined naturally 
in the way we have described for MPS. 

Let us now consider numerical examples for different tar- 
get states 10) to demonstrate the feasibility and efficiency 
of the proposed algorithm. We start with ground states of 
nearest-neighbor Hamiltonians on a chain, i.e., the |0) = \gs) 
are completely determined by all the reductions to two ad- 
jacent spins as outlined above and the above algorithm not 
only produces states that match the reduced density matrices 
of the ground states but, in fact, states that are themselves 
close to the ground states. Among ground states of one- 
dimensional nearest-neighbor Hamiltonians the critical ones 
are the most challenging to approximate by MPS as they vio- 
late the entanglement-area law |23 | and we test our algorithm 
for such an example: the critical Ising model. In order for the 
local operators in the Hamiltonian not to be exactly the ones 
that are measured, we also rotate the Ising model locally by 
7r/4 around each spin axis. This model is solvable and in order 
to show that we do not consider a pathological case, we also 
consider one-dimensional random Hamiltonians of the form 

i=l 

(i) (i) 

where the f- \ f-^^ act on spin i and z + 1, respectively, and 
are hermitian matrices with entries that have real and imagi- 
nary part picked from a uniform distribution over [—1,1]. For 
each Hamiltonian, we first determine the ground state \gs) ex- 
actly (i.e., the target state 10)) and its reductions and then com- 
puted the fidelity \{gs\yn)\'^ after n iterations of the MPS-SVT 
algorithm, see Fig.[T] 

Our method is of interest for all situations in which stan- 
dard tomography will not be feasible. This is the case for the 
verification of state preparation in experiments with too many 
particles. An example is the recent ion trap experiment |4 1 for 
the preparation of W-states, |0) = (|10 • • • 0) + |010 • • • 0) + 
• • • + |0 • • • 01))/>/]V, that were limited to 8 qubits princi- 
pally because the classical postprocessing of data became pro- 
hibitive for longer chains. Here we demonstrate the efficiency 
of our approach (we are not limited to few ions and demon- 
strate convergence for up to 20 ions - even higher number of 
ions are easily accessible due to the MPS alteration of the SVT 
method) by illustrating how one would postprocess experi- 
mentally obtained reduced density matrices to guarantee the 
generation of 1 0) or a state very close to it. We mimic experi- 
mental noise by adding Gaussian distributed random numbers 
with zero mean to the pk. After initializing the MPS algo- 
rithm with the MPS representation of 1 0) and Yq = R, we use 



3 




FIG. 1: Fidelity fN,n = |((/)|2/n)|^ as a function of the number of spins N and iterations n of the MPS-SVT algorithm for different target 
states 10). Left: Ground state of the locally rotated critical Ising model, fN,n (surface, left axis) and 1/(1 — fN,n) (lines, right axis), showing 
that, for fixed system size, 1 — fN,n decreases as Middle: i/l — fN,n as a function of N for the ground state of the locally rotated 

critical Ising model (lines, right axis, from top to bottom n = 5, 10, 50, 100, 200, 500, dots are obtained by exact numerical diagonalization of 
the Hamiltonian and the Yn) and random Hamiltonians (1000 for each N) as described in the main text after n = 5 iterations (densities, left 
axis, arrows indicate the mean). In both cases the scaling for fixed n of 1 — fN,n is better than ^ N"^. Right: W state as described in the text 
for 4000 MPS-SVT iterations. Plot shows | ((/)|2/n) |^ as a function of the number of ions, for no noise (dots) and Gaussian noise (densities 
obtained from 100 realizations for each N, arrows indicate mean) with a standard deviation of 0.005 (even N) and 0.01 (odd N). 



Y.^\pk - {yn\Pk\yn) \ ^ figure of merit for conver- 
gence, i.e., after a given number of iterations, we pick the \yn) 
with minimum Xn • The result of such a procedure is shown in 

So far we have presented the method for pure states and 
one-dimensional systems. The SVT algorithm as described 
above works for higher-dimensional systems as well and may 
be made efficient by adopting higher dimensional analogues 
of MPS methods as outlined above. The extension to mixed 
states is also straightforward as its treatment can be reduced to 
that of pure states by using the fact that every mixed state on N 
qubits can be purified to a pure state on 27V qubits. Hence we 
may ask for a globally pure state on 2N qubits that matches 
the reduced density on all contiguous sites of k qubits on the 
first N qubits. While the reduced density matrices do not 
uniquely determine the mixed state, approximations of bet- 
ter and better quality can be obtained by increasing k. As an 
example, suppose the state is the Gibbs state corresponding to 
a A:-local Hamiltonian H, i.e., the state q minimizing the free 
energy 

tr[^^]-T^(^). (4) 

The first term is, as before, for a /c-local Hamiltonian, deter- 
mined by the reduced density matrices. The entropy of the to- 
tal state however can only be learnt exactly from the complete 
density matrix. However, for essentially all reasonable physi- 
cal systems, the entropy density lim/e^oo S{trk+i^...{p)) /k in 
the thermal state of a Hamiltonian exists ll24ll and as a con- 
sequence the total entropy of the state can be estimated effi- 
ciently from the knowledge of reduced density matrices. 



Our algorithm described above may also be adapted 
straightforwardly to determine hypothesis states in recently 
proposed algorithms for quantum learning l|25ll . Here, a small 
given set of randomly chosen observables is measured and 
on the basis of the measurement outcomes a quantum state 
closely approximating the measured expectation values needs 
to be found. This state will, with large probability, predict the 
expectation values of all observables. Present approaches to 
determine such states are based on semi-definite programming 
and are therefore inherently non-scalable, they are limited to 
perhaps 12 qubits 1251 . 

Our algorithm will be essential for efficient tomography and 
verification of medium to large scale quantum information de- 
vices. Already today they are beginning to reach scales for 
which standard tomography is not feasible anymore. Further- 
more, it may also be applied to problems in condensed matter 
physics where the system has too many components to achieve 
tomography by standard means. 

Hence our combination of singular value thresholding with 
the matrix product state representation is expected to become 
a useful tool in a wide variety of physical settings and algo- 
rithms. 
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Appendix 

In singular value thresholding IH [51 one seeks a solution 
for the minimization of the trace norm of a matrix X subject 
to some linear constraints, i.e. 

minimize tr|X| 
subjectto7^Q(X) =7^q(M) 



where Vq{M) is a matrix whose elements are non-zero only 
on the entries belonging to the index set For a value r > 
and a sequence {Sk}k>i one inductively defines 

X(^) =shrink(r(^-^\r) 

where, in standard SVT, 

shrink(y, r) = [/diag({max{0, - r}})V'^ (5) 

with the singular value decomposition Y = Udmg{{ai})V^ . 

If it is our goal to reconstruct pure states compatible with 
given reduced density matrices we may adapt SVT and in 
the process make it suitable for application to matrix product 
states. To this end we introduce a small but crucial variation 
of the shrink operation. Rather than introducing a threshold 
r we retain only the largest singular value ai and the corre- 
sponding matrix U diag(max^ • • • 0) If the target state 
is pure, i.e., a matrix with rank equal to one, this can be ex- 
pected to converge rapidly, an expectation that is confirmed 
by extensive numerics. 

While Vn{M — X*^^^) is not itself positive, initializing the 
recursion with a positive operator, e.g. a pure state, and choos- 
ing sufficiently small 6^ will ensure that in the second step of 
the recursion a matrix is generated whose largest eigenvalue 
is positive and equal to the largest singular value. 

Crucially, this largest eigenvalue and corresponding eigen- 
state can then be computed efficiently via the maximization 
of the expectation value of a matrix product state solving 
the multi-quadratic optimization problem by a succession of 
quadratic optimization problems each of which can be solved 
via a generalized eigenvalue problem IT61 . Furthermore, 
Vn{M -X'^^^) in the recursion is a sum of a linear number 
of Paulistrings whose expectation values in a matrix product 
state may be obtained efficiently. 

Reduced density matrices are obtained by measuring the 
set of strings of Pauli operators Pk, k = 1, . . . , iC. Then 
T^n{Y) = Ef=i tr[r A] and is hermitean. 



